Abstract 


A multiple scales approach is used to approximate the effects of nonparallelism and 
streamwise curvature on the stability of three-dimensional disturbances in incompressible 
flow. The multiple scales approach is implemented with the full second-order system of 
equations. A detailed exposition of the source of all terms is provided. 


1 Introduction 

Multiple scales analysis is a useful tool in fluid mechanics. It is often used to derive weakly 
nonlinear and weakly nonparallel corrections to eigenvalue problems that arise in the study of 
the stability of laminar flows. Although the general procedures are well-known [1, 2], the novice 
can easily become bogged down in the derivation of the appropriate equations. This note is 
intended to clarify the detailed steps that are needed in the analysis. 

Previous analyses of the weakly nonparallel and curvature effects have been based on a 
reduction of the governing equations to a first-order system of equations [1, 2, 3]. This reduction 
in order was accomplished with the introduction of many additional variables; although these 
additional variables are mathematically well-defined, the physical importance of their adjoint 
counterparts, which are required for the solution of the problem, is nonintuitive. Historically, 
the first-order formulation of stability equations was popularized at a time when the numerical 
methods for the solution of linear eigenvalue problems were first being developed. With the 
addition of orthonormalization, previously developed marching techniques for first-order systems 
of differential equations could be applied to eigenvalue problems that were formulated as first- 
order systems. As computing power has increased, matrix methods, which are more suitable for 
obtaining the global eigenvalue spectrum, have largely superseded marching methods. However, 
to compute nonparallel effects, researchers have continued to use the first-order formulation 
[4, 5]. We show that the multiple scales analysis can also be performed easily in the context of 
the original equations, with the original, physically meaningful variables. 

Masad and Malik [4] previously used multiple scales analysis to include surface curvature 
for the case of flow over a swept circular cylinder, hence a cylindrical coordinate system was 
appropriate for their application. Our goal is to extend the analysis to a wider range of bound- 
ary layers, hence the analysis is presented in surface-fitted coordinates. The mathematics are 
presented in section 2; numerical verification is provided in section 3. 


2 Mathematical Model 

For expository purposes, we focus on the linear stability of a three-dimensional, incompressible 
boundary layer. However, apart from the specification of the operators involved, most of the 
analysis is directly applicable to compressible flows as well. All lengths are nondimensionalized 
by a length scale L*; all velocities, by a velocity scale U *; time, by the ratio L* jU*\ and pressure, 
by p*U* 2 , where p* is the density. The flow Reynolds number is R = U*L* /v * , where v* is the 
dimensional kinematic viscosity. The streamwise direction is denoted by x] the wall-normal 
direction, by y\ and the spanwise direction, by z. The base flow 

Qo = (U 0 ,eV 1 ,W 0 ) (1) 
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consists of 0(1) streamwise and spanwise components and a small 0(e) surface-normal com- 
ponent. For the zero-pressure-gradient boundary layer, only the surface-normal component of 
mean velocity has an 0(e) term; other velocity component corrections enter at higher orders. 
To simplify the exposition, we assume that the mean flows considered here can be approximated 
in a similar manner. All mean-flow components are invariant in both time t and spanwise di- 
rection. The mean flow varies rapidly in the surface-normal direction and changes slowly in the 
streamwise direction. The slow variation in the streamwise direction suggests that a slow scale 
x\ = ex should be introduced. We consider the disturbance quantities 


q (x,y,z,t) = [A(zi)qo(a:i, 2 /) + eqi(aq ,y) + ...}e l6 


(2) 


where q = [u,v,w,p] T represents the vector of the streamwise, surface-normal, and spanwise 
disturbance velocities, and the disturbance pressure, respectively. The variable A(x\) represents 
the amplitude of the disturbance, and qo (%i,y) denotes the vector of the normalized quasi- 
parallel eigenfunctions. The 0(e) term in Eq. (2) denotes the nonparallel correction to the 
disturbance. The argument 0 of the exponential function is defined in terms of the streamwise 
and spanwise wave numbers and the circular frequency of the disturbance. That is, 
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Logarithmic differentiation of Eq. (2) with respect to x yields the growth rate of the disturbance: 


— Imag (a) = — Imag [ao(aq) + ea i] + 0(e 2 ) 


( 4 ) 


The term —Imag [ao(aq)] is the quasi-parallel growth rate. The ea.\ term is given by 

_ _ . f 1 dA(x i) 1 dq 0 (x 1 ,y 1 y 

€ai 1 [A dx t q 0 (x 1 ,y 1 ) dx t 


( 5 ) 


of which the negative of the imaginary part represents the nonparallel growth of the amplitude 
function A and the change in shape of the quasi-parallel eigenfunction. Unlike the quasi-parallel 
case, the growth rate in the nonparallel context depends on both the particular disturbance 
quantity qo and the surface-normal location y\. Here, we choose qo{x\,yi) to represent uq (i.e. , 
the streamwise velocity component of qo at the surface-normal location at which uq is at a 
maximum). This location varies with x\. 

To obtain the nonparallel equations, the mean and disturbance quantities are substituted 
into the Navier-Stokes equations. Equations for the mean flow Qo ( x i,y) are solved separately. 
The disturbance equations are obtained by hnearizing and equating like powers of e. The 0(1) 
set of equations can be written as a multivariable second-order form of the Orr-Sommerfeld 
equation: 

L os q 0 = 0 (6) 

with boundary conditions, 

uq = v 0 = wq = 0 at y = 0 (7) 


and 


M o, v o, w o,Po 0 as y oo 


(8) 
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With the notation D = d/dy , k 2 = (oq + /3 2 ), and 5 = k 2 + iR(aolIo — u) — D 2 , the Orr- 
Sommerfeld operator for incompressible flow is 


/ S RDUo 0 ia 0 R \ 

0 S 0 RD 

0 0 5 -i(3R 

\iao D if3 0 / 


(9) 


Note that L os is a matrix operator that can be expressed as 

L os = Lq + (iao)Li + (ioo) 2 L2 (10) 


where Lo, Li, and L 2 are independent of op and depend only on the mean flow and the stability 
parameters u and [3. The operator Lo can be thought of as arising from those terms in the 
primitive- variable equations that do not contain any x derivatives. Similarly, the operators Li 
and L 2 are the coefficient matrices for terms that contain only first and second x derivatives 
respectively. The multiple scales analysis indicates that the first derivative of the assumed 
zeroth-order solution Aq^e 16 with respect to x can be written as 


d(Aq 0 e l 

dx 


( dA <9q 0 

ia 0 Aq 0 + e - — q 0 + A - — 

\OX\ OX\ 


( 11 ) 


and the second derivative with respect to x can be written as 

d 2 (Aq 0 e l 


dx 2 


= (ia 0 y Aq 0 e l 


+ e 


(9A f 


2ia ° i + + 1 

\ux 1 ox \) Ox 1 J 


e %d + 0 (e 2 ) 


( 12 ) 


The 0(e) problem is obtained by substituting Eqs. (2), (11), and (12) into the Navier-Stokes 
and continuity equations. Thus we find, 


-r 9A^ t , O , <9qo dia 0 ^ T 

L os qi = — -v — M 0 qo — A I M 0 -v h — Miq 0 + Nq 0 

OX\ \ OX\ OX\ 

where Mq = Li + 2iapL2; Mi = L 2 ; and N contains nonparallel mean-flow terms: 


(13) 


N = 


/-^(f f + ViD) 0 

0 -RiDVx+VxD) 

-R^ 0 

da 7i 

Vo 0 


°\ 

0 


-RV-iD 0 

0 0 J 


(14) 


The boundary conditions on the 0(e) problem are 

u\ = v\ = w\ = 0 at y = 0 

and 

U\ , Vi , vji , pi — ^ 0 as y — ► 00 


(15) 


(16) 


Because the homogeneous part of the O(e) problem is the same as the 0(1) problem, the O(e) 
problem has a solution if and only if a special condition (known as the solvability condition) is 
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satisfied [6]. For the solvability condition to be satisfied, the right-hand side of Eq. (13) must be 
orthogonal to every member within the null space of the adjoint of L os . Hence, if we dot-multiply 
the adjoint solution with the right-hand side of Eq. (13), integrate over the domain, and set 
the result equal to zero, we obtain the expression 


OO / \ 

/ • (Mogfc + t^M iqo + Nq 0 ) dy 

OO 

/ • (M 0 q 0 )dy 


(17) 


The imposition of the solvability condition ensures that no secular terms exist in the solution. 
(For a discussion of secular terms, see, for example, Ref. [7].) 

Typically, the adjoint problem is obtained by first transforming Eq. (6) into a set of first- 
order differential equations [1, 2]. However, here we show that the adjoint problem is easily 
obtained from the governing equations in their natural form. First, premultiply Eq. (6) by q^ 
and then integrate over the surface-normal direction. The adjoint operator is obtained by using 
integration by parts and the linear- algebra identity, (yPx) T = x T P T y T , where the superscript 
T denotes transpose. The adjoint operator can be written as 

= A t D 2 - B t D + (C - DB t ) (18) 


where 

L os = AD 2 + BD + C (19) 

The component matrices A, B, and C are defined as 
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/ k 2 + iR(a 0 Uo — w) 
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C = 


0 


k 2 + iR(a 0 Uo — w) 


( 20 ) 


ia 0 R \ 
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k 2 + iR(aolIo — u) —i(3R 

0 


) 


0 0 

iao 0 i(3 

Note that A = A T ; and for incompressible how, D B T = 0 so that the adjoint operator can be 
simplified. 


Two terms on the right-hand side of Eq. (13) must be determined before the solvability 
condition can be applied. These terms are the x\ derivatives of iao an d qo- Note that the 
partial derivative of the 0(1) equation Eq. (6) with respect to x\ must satisfy 


$L os qo 9Tj 0S d Qo d L os d{iao ) T 9qo „ 

O Ti QO ' o/ • \ o Qo "T -L'os o U 

OX\ ^^0 OX\ UylQLQ) OX\ OX\ 


( 21 ) 
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After the known terms are moved to the right-hand side, Eq. (21) can be rewritten as 


L 


OS 


dq 0 

dx\ 


' dL 0S (9Q 0 dL 0S d(ia 0 ) 

. <9Qo dx 1 d(ia 0 ) dx 1 


( 22 ) 


As in Eq. (13), this expression has a nontrivial solution if and only if the right-hand side is 
orthogonal to every solution of the homogeneous adjoint problem. The solvability condition 
that determines d(iao)/ dx\ is obtained in the same manner as described previously to obtain 
the result in Eq. (17). With d(iao)/dxi being known, equation (22) can then be solved for 
dqo/dxi. Finally, both terms are available to evaluate using Eq. (17). 

To obtain the complete 0(e) correction to the growth rate in Eq. (5), one must make a 
specific choice for dqo/dxi, i.e. , chose the physical quantity qo of interest and the surface normal 
location, y\ where the growth rate is to be determined. The term dqo/dxi accounts for the 
slow change in the shape of the eigenfunction with downstream distance. Two commonly used 
choices for qo(x\, y\) are (1) uq at the y\ value at which uq is at its maximum (referred to as the 
K mas choice below) and (2) uq at a predefined y\ location. In both cases, dqo/dxi from Eq. (22) 
defines the value of dqo/dxi after the selection of variable and surface-normal location has been 
made. 

As shown by Masad and Malik [4], weak surface curvature can be treated as a small per- 
turbation to the problem without curvature. The modifications that are needed to incorporate 
curvature effects into the equations above are minor and are easily implemented. First, the 
curvature introduces terms that are proportional to ^-j^- where r is the local radius of curvature. 
For a large r, these additional terms are simply proportional to the smaU surface curvature 
ck = 1. Because the terms are of 0(e), they can be easily added to the matrix N in Eq. (13). 
In body-fitted coordinates, the necessary additional terms can be surmised from Lin and Reed 
[8] and are given as matrix N c : 


/ RV\k RUok 0 0 \ 

-2 RU 0 k 0 0 0 

0 0 0 0 

\ 0 k 0 0 / 


(23) 


Secondly, the curvature modifies the x derivative terms so that they are proportional to 
For large r, can be approximated by 1 — = 1 — eny. If we keep only those terms that are 

first order in ck, then the additional O(e) terms are M c qo, where 


M c = ia 0 ny Li 


(24) 


Hence, the modified O(e) equation that must satisfy the solvability condition is 


dA 

L os qi = -x — M 0 q 0 - A 

OX i 


M o |qo + ^ Miqo + (M c + N + N c ) qo 


dx-i 


dxi 


(25) 


This equation is solved with the same procedures that are required for Eq. (13). 

The numerical solution of these equations involves the discretization of the equations with 
the use of Chebyshev polynomials. A staggered grid is used for the pressure variable and, hence, 
the continuity equation. An iterative procedure is used to determine the quasi-parallel eigenvalue 
and the eigenfunction. Direct solves are used to solve the systems of equations. The computer 
code used for the calculations is a modified version of SPECLS [9]. 
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3 Verification 


Figure 1 shows pointwise comparisons of the nonparallel results with the multiple scales results 
obtained by El-Hady [2] for incompressible flow at R = 1000, for a variety of nondimensional 
frequencies (F = 2ir f*v* /U* 2 , where /* is the dimensional frequency). The surface-normal axis 
represents the change in the growth rate (— Imag(eai)) due to the nonparallel effects at the 
wave number of maximum quasi-parallel growth rate. Both our results and those of El-Hady 
are based on the streamwise component of velocity at the surface-normal location at which the 
linearized solution is maximum. The only significant difference between the approaches used 
involves the order of the system of equations. El-Hady used a first-order system of equations; 
we use the second-order system. The agreement of the two sets of results strongly suggests that 
both approaches are correctly implemented. 

A similar comparison was made to test the implementation of the weak curvature terms. 
In this case, we consider the incompressible flow over a circular cylinder that is swept at an 
angle of 60.5° relative to the free stream. The Reynolds number, based on the free-stream speed 
and the radius of the cylinder, is 2.0 x 10 6 . We consider only steady ( F = 0) disturbances 
with a spanwise wave number of 0.8 (normalized with the radius). Weakly nonparallel effects 
are neglected. The comparisons are made with earlier data from Masad and Malik [4]. Masad 
and Malik [4] also used a multiple scales analysis to incorporate the curvature effect; however, 
because their work focused only on the cylinder case, they used a cylindrical coordinate system 
to derive the perturbation equations. They also used a somewhat less accurate numerical scheme 
that employed finite differences. However, as shown below, the results that we obtained are in 
excellent agreement with their results. In Fig. 2(a), the real part of the streamwise wave number 
a is shown both with and without curvature for various angles 0* relative to the stagnation 
point on the cylinder. The effect of the curvature terms is relatively small here, but the results 
obtained with our multiple scales analysis essentially duplicate those obtained by Masad and 
Malik [4]. A comparison of the results for the imaginary part of a in Fig. 2(b) show the same 
excellent agreement. 
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Figure 1. Nonparallel correction to growth rate at R = 1000 in a Blasius boundary layer: 
comparison with El-Hady [2]. 
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(a) Real part of a. 



(b) Imaginary part of a. 


Figure 2. Multiple scales approach for weak surface curvature: comparison with data from 
Masad [4]. R = 2.0 - 10 6 ; flow angle 60°; F = 0.0. o current quasi-parallel; □ current with 
curvature; A Masad quasi-parallel; o Masad with curvature. 
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